Transcriptome analysis reveals differential transcription in tomato (Solanum lycopersicum) following inoculation with Ralstonia solanacearum

Tomato (Solanum lycopersicum L.) is a major Solanaceae crop worldwide and is vulnerable to bacterial wilt (BW) caused by Ralstonia solanacearum during the production process. BW has become a growing concern that could enormously deplete the tomato yield from 50 to 100% and decrease the quality. Research on the molecular mechanism of tomato regulating BW resistance is still limited. In this study, two tomato inbred lines (Hm 2–2, resistant to BW; and BY 1–2, susceptible to BW) were used to explore the molecular mechanism of tomato in response to R. solanacearum infection by RNA-sequencing (RNA-seq) technology. We identified 1923 differentially expressed genes (DEGs) between Hm 2–2 and BY 1–2 after R. solanacearum inoculation. Among these DEGs, 828 were up-regulated while 1095 were down-regulated in R-3dpi (Hm 2–2 at 3 days post-inoculation with R. solanacearum) vs. R-mock (mock-inoculated Hm 2–2); 1087 and 2187 were up- and down-regulated, respectively, in S-3dpi (BY 1–2 at 3 days post-inoculation with R. solanacearum) vs. S-mock (mock-inoculated BY 1–2). Moreover, Gene Ontology (GO) enrichment analysis revealed that the largest amount of DEGs were annotated with the Biological Process terms, followed by Cellular Component and Molecular Function terms. A total of 114, 124, 85, and 89 regulated (or altered) pathways were identified in R-3dpi vs. R-mock, S-3dpi vs. S-mock, R-mock vs. S-mock, and R-3dpi vs. S-3dpi comparisons, respectively, by Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway analysis. These clarified the molecular function and resistance pathways of DEGs. Furthermore, quantitative RT-PCR (qRT-PCR) analysis confirmed the expression patterns of eight randomly selected DEGs, which suggested that the RNA-seq results were reliable. Subsequently, in order to further verify the reliability of the transcriptome data and the accuracy of qRT-PCR results, WRKY75, one of the eight DEGs was silenced by virus-induced gene silencing (VIGS) and the defense response of plants to R. solanacearum infection was analyzed. In conclusion, the findings of this study provide profound insight into the potential mechanism of tomato in response to R. solanacearum infection, which lays an important foundation for future studies on BW.

or honeycomb. When the diseased plant is cut open, the fibrous tube tissue inside will turn brown. When the transverse cut of the diseased plant base is pressed hard, the yellow-white bacterial mucus will flow out of the fracture, which is known as "bacteria pus" 12,13 . Bacterial wilt, known as plant cancer, seriously affects the yield and quality of crops. It has strong variability and soil-borne characteristics. The current traditional control methods, such as breeding resistant varieties, crop rotation and chemical control have some limitations [14][15][16] . Therefore, it is necessary to have a comprehensive and detailed understanding of plant-pathogen interactions during the progression of BW 17 .
The current research on tomato resistance to BW mainly focuses on the genetic mechanism of resistance [18][19][20][21] ; identification 20,22,23 and screening 13,[24][25][26][27][28][29] of molecular markers related to the disease resistance genes and other aspects. There are only a few studies available on the gene regulation of BW resistance [30][31][32][33][34][35][36][37] , and the molecular mechanism of tomato resistance-related genes regulating BW remains unclear. Thus, it necessitates performing the research on the response of tomato plant to BW. In recent years, with the application of transcriptome sequencing technology, several genes and miRNA functions have been identified. Transcriptome sequencing has been widely used in basic research, molecular breeding, pathogen-host interaction mechanism, comparison of resistance genes between susceptible and disease-resistant varieties, and biocontrol factors inducing plant disease-resistance mechanisms [38][39][40] . Transcriptomic technique has been used in many studies to identify the molecular mechanism of R. solanacearum resistance in a plethora of plant species, including Solanum dulcamara 41 , Arachis hypogaea 42 , Solanum commersonii 43 , Solanum melongena 44 , Capsicum annuum 45 , Solanum tuberosum 46 , and Nicotiana tabacum 47 . In tomato, the in-depth transcriptome data is available for its interaction with R. solanacearum 15,17,[48][49][50] . French et al. 48 analyzed the genome-wide transcriptional response of roots of resistant and susceptible tomato plants at multiple time points after inoculation with R. solanacearum and identified the molecular basis of this resistance. Furthermore, they determined the role of auxin signaling and transport pathways in resistance against R. solanacearum by functional analysis of an auxin transport tomato mutant. Jiang et al. 49 investigated the root transcriptome profiles between silicon (Si)-treated (+ Si) and untreated (− Si) tomato plants at different days post-inoculation with R. solanacearum by using RNA-seq technology. They also determined the content of hormones including salicylic acid (SA), jasmonic acid (JA), and ethylene (ET). Finally, they suggested that Si enhanced BW resistance of tomato via several signaling pathol.ways. The molecular mechanism of tomato resistance-related genes regulating bacterial wilt resistance is still unclear. Therefore, it is necessary to screen out the genes related to tomato BW resistance through transcriptome analysis.
In this study, we used RNA-seq technology to analyze the transcriptome of two tomato inbred lines resistant to BW Hm 2-2 (R) and susceptible to BW BY 1-2 (S) before and after inoculation with R. solanacearum. The two tomato inbred lines are special and the inductive properties are stable. In addition, studies comfirmed that Ralstonia solanacearum can multiply in plant stems. In this study, tomato stems were collected and studied, while tomato roots were studied in most previous studies. This study aimed to determine the molecular mechanism involved in the tomato response towards R. solanacearum, and provide a theoretical foundation for future research in BW.

Results
Phenotypic characterization after inoculation with R. solanacearum. At 3 dpi (days post-inoculation), plants exhibited different phenotypic symptoms. As shown in Fig. 1, leaves of the resistant plants (Hm 2-2) and susceptible (BY 1-2) plants inoculated with sterile water showed no obvious symptoms. However, one  Analysis of RNA-seq data. To obtain a global overview of the transcriptome relevant to BW stress conditions in the resistant and susceptible tomato plants, the cDNA libraries from stem samples of Hm 2-2 and BY 1-2 plants with mock-inoculation and pathogen-inoculation were sequenced on the Illumina HiSeq2500 platform, separately. After removal of the adaptors and filtration of low-quality reads, an average of 3.3 × 10 7 clean reads were obtained for each sample, and clean reads from all twelve libraries were mapped to the tomato genome ITAG3.2. The coverage of overall mapped reads to all the samples ranged between 96-97% and the percentage of uniquely mapped reads ranged between 76 and 85%. Moreover, the percentage of bases with Q20 (high sequencing quality) was close to 98%. Detailed sequencing and mapping statistics are given in Supplementary Table S1. As shown in Supplementary Fig. S1, gene coverage ranged from 80 to 100%, accounting for approximately 75% of the total genes ( Supplementary Fig. S1a). An experimental repeatability test was performed, the results of which indicated that the results between replicates were comparable ( Supplementary Fig. S1b). These results showed that the transcriptome sequencing quality was appropriate for further analyses.
Identification of DEGs. Differential expression analysis was performed by the DESeq2 package between two groups (treatments and control). The genes with FDR < 0.05 and |log2FC|≥ 1 were considered DEGs. The number of DEGs was different in each of the four comparisons (Fig. 2a). The most DEGs were found in S-3dpi vs. S-mock (1087 and 2187 significantly up-and down-regulated genes, respectively), followed by R-3dpi vs. R-mock (828 and 1095 significantly up-and down-regulated genes, respectively) and R-mock vs. S-mock (386 and 258 significantly up-and down-regulated genes, respectively), and the least DEGs were seen in R-3dpi vs. S-3dpi (337 and 263 significantly up-and down-regulated genes, respectively). In R-3dpi vs. R-mock and S-3dpi vs. S-mock, the total number of down-regulated genes was greater than that of up-regulated genes. While in R-mock vs. S-mock and R-3dpi vs. S-3dpi, the total number of up-regulated genes was greater than that of down-regulated genes. The volcano plots showed that the number of up-and down-regulated genes had a clear distribution pattern in all four comparisons (Fig. 2b). Moreover, the distribution pattern of down-regulated genes in S-3dpi vs. S-mock was much higher than that in R-mock vs. S-mock and R-3dpi vs. S-3dpi. Although the numbers of up-and down-regulated genes in R-mock vs. S-mock and R-3dpi vs. S-3dpi were lower compared with those of the other two comparisons, the distribution patterns were very similar. These results clearly showed the global gene expression patterns between different comparisons. The heat map analysis based on hierarchical clustering analysis of DEGs in four comparisons showed that DEGs in different comparisons showed different expression trends. Among the top 50 highly correlated genes across samples, 31/19 clustered together to reveal higher up-/down-regulating trends in R-3dpi vs. R-mock and S-3dpi vs. S-mock comparisons. Most of the top 50 highly correlated genes exhibited little difference between R-mock vs. S-mock and R-3dpi vs. S-3dpi (Fig. 3). The expression of the top 50 DEGs in different comparisons is listed in Supplementary Table S3. SNP and InDel annotations. Variation detection based on SNPs and InDels of transcriptome sequencing was performed by the software GATKv3.4-46 (Fig. 4). Two types of nonsynonymous single nucleotide variation (SNV) and synonymous SNV represented dominant functional variations among the nine types of functional variations (Fig. 4a). The remaining types were frameshift deletion, stop gain, unknown, non-frameshift deletion stoploss, frameshift insertion, and non-frameshift insertion. Furthermore, these SNPs were distributed in different locations, with dominant distribution in intronic regions, followed by intergenic and exonic regions (Fig. 4b). We also identified all the possible mutations in this study. The two dominant types were transition and transversion (Fig. 4c), where transition accounted for 61.44%, while transversion accounted for 38.56%. These results suggest systematic, comprehensive, transcriptional regulation in tomato that orchestrates the response to BW.
GO enrichment analysis of DEGs. GO enrichment analyses were performed for all significant DEGs (Fig. 5). Different comparisons showed similar distribution patterns in terms of the number and type of enriched GO terms, which could be divided into three main functional groups, including 22 biological processes, 14 molecular functions, and 15 cellular components (Fig. 5a). The largest amount of DEGs were annotated with the biological process terms, with the cellular process, metabolic process, and single-organism process accounting for the largest number of DEGs. In molecular function terms, the majority of DEGs were associated with binding and catalytic activities. Other significantly abundant cellular component terms were cell, cell part, membrane, organelle, and membrane parts (Supplementary Table S4). However, the enrichment level (Q-value) for each functional comparison varied (Fig. 5b). Notably, other functional comparisons that were significantly enriched were involved in different cellular component pathways, such as extracellular region, cytoskeleton, microtubule cytoskeleton, external encapsulating structure, and cell periphery processes (Fig. 5b).

Plant-pathogen interaction pathways.
Previous studies have concluded that increased Ca 2+ -dependent cyclic nucleotide-gated channels (CNGCs) play a key role in plant response to pathogens and pathogen-associated molecular pattern (PAMP) signals 51 . As shown in Supplementary Fig. S2a, the expression levels of CNGCs and HSP90 were up-regulated in R-mock vs. S-mock, but after R. solanacearum inoculation, their expression levels were significantly down-regulated ( Supplementary Fig. S2b) in R-3dpi vs. S-3dpi. Furthermore, the CERK, MEKK1, Rboh, and PBS1 genes were down-regulated, but MYC2, PR1, and RIN4 were up-regulated in R-3dpi vs. S-3dpi. These results can help understand the resistance mechanism of tomato to R. solanacearum.
Plant hormone signal transduction pathways. Most of the plant hormone signal transduction pathways were represented by DEGs in both genotypes with only a few differences. Auxin responsive AUX/IAA and SAUR were found up-regulated in both R-mock vs. S-mock and R-3dpi vs. S-3dpi. Cytokine responsive B-ARR was found up-regulated in R-mock vs. S-mock but not significantly changed in R-3dpi vs. S-3dpi comparison. Moreover, AHP was found down-regulated in R-3dpi vs. S-3dpi but no significant change in AHP expression was found in R-mock vs. S-mock. The TGA transcription factor, a regulator of NPR1, was found up-regulated  Fig. S3).  Fig. S4). EBF1/2 genes in ethylene-responsive defense response were found down-regulated in both R-mock vs. S-mock and R-3dpi vs. S-3dpi comparisons, but interestingly, ChiB that was up-regulated in R-mock vs. S-mock showed down-regulation in R-3dpi vs. S-3dpi. Moreover, EIN3/EIL was found up-regulated in R-mock vs. S-mock but not significantly changed in R-3dpi vs. S-3dpi. Few plant-defense-related components, such as PR1 was found up-regulated, whereas EBF1/2 were commonly down-regulated in both R-mock vs. S-mock and R-3dpi vs. S-3dpi. Similarly, FLS2 was found down-regulated in both R-mock vs. S-mock and R-3dpi vs. S-3dpi. MPKKK17/18 were found up-regulated in R-3dpi vs. S-3dpi but not significantly changed in R-mock vs. R-mock. Furthermore, cell death-related RbohD was found downregulated in R-3dpi vs. S-3dpi but not significantly changed in R-mock vs. S-mock ( Supplementary Fig. S4).

Validation of RNA-seq data by qRT-PCR analysis.
To verify the reliability of Illumina sequencing results, eight genes were randomly selected for qRT-PCR analysis with three biological replicates. The expression levels of these genes were normalized to the constitutively expressed Actin gene.  Fig. 7, and the primers for qRT-PCR are given in Supplementary  Table S2. Overall, the results obtained from qRT-PCR and RNA-seq showed the same up-and down-regulation trends, which suggested that the RNA-seq results were reliable.   Fig. 8a. Figure 8b shows the expression pattern of WRKY75 in VIGS plants. We evaluated the pathogenicity of the WRKY75-silencing plants in comparison with the wild-type Hm 2-2 plants (Fig. 8c). The result showed that the disease rating of WRKY75-silencing plants was greatly increased compared with the wild-type plants.

Discussion
RNA-seq technology, namely "whole transcriptome shotgun sequencing" technology, is an effective method for comprehensive analysis of the entire transcriptome through high-throughput sequencing 52 . Compared with the conventional hybridization technique, RNA-seq technology, without advanced probe design, can test all transcription information for any species at the single nucleotide level, and it has gradually become a new technology of genome and transcriptome research due to its advantages including high-throughput, easy to operate, can be quantitative, and low cost. It is widely used in differential gene detection, new gene identification, and gene   56 performed RNA-seq analysis between two genotypes containing PMR-1 (resistant) and Pusa Vishal (susceptible) mungbeans in response to yellow mosaic virus. The resistance to mungbean yellow mosaic India virus (MYMIV) showed a very complicated gene network, starting from the production of general PAMPs, and activating various signaling cascades such as brassinosteroid (BR), jasmonic acid (JA), and kinases. Finally, the expression of specific genes (such as R-gene proteins, virus resistance, and PR-proteins) leads to disease resistance response. In this study, we identified 1923 and 3274 DEGs in R-3dpi vs. R-mock and S-3dpi vs. S-mock comparisons, respectively, by using RNA-seq technology. The obtained results revealed that these DEGs were involved in plant metabolism, signal transduction, synthesis of secondary metabolites, genetic information processing, responses to environmental stimuli, self-immunity, and other life activities. These results showed that the quality of the transcriptome sequencing output and assembly that we obtained meets the requirements of transcriptome analysis. The use of RNA-seq www.nature.com/scientificreports/ technology can yield very comprehensive and abundant transcription information, and that can provide us with abundant resources and paths for future research. Plant hormones are a class of organic substances produced by the plant itself. They play an important role in regulating various crucial activities including growth, development, senescence, and death at very low concentrations. Relevant studies have shown that SA, JA, and ET are the crucial signal molecules that induce plants to produce defense responses [57][58][59] . In addition, BRs, auxins, and cytokinins (CTKs) are also taking part in the defense response of plants to pathogens [60][61][62] . As a gaseous hormone, ET not only has a wide-ranging regulatory effect on plant growth and development but also participates in plant disease resistance and defense response [63][64][65][66] . Tezuka et al. 67 reported that the rice ethylene response factor (ERF) OsERF83 positively regulates blast resistance by regulating the defense-related genes' expression in rice. In this study, the expression of EBF1/2 was up-regulated in the S-3dpi vs. S-mock comparison (Supplementary Fig. S3b). SA is a type of phenolic hormone widely distributed in plants. It participates in various physiological processes of plants and plays a key role in forming defense responses against pathogenic bacteria [68][69][70] . It has been proved that after plants are infected by pathogenic bacteria, the SA that increases exponentially in the body can induce plants to produce hypersensitive response (HR) and systemic acquired resistance (SAR) so that plants exhibit disease resistance 71 . The SA-mediated disease resistance signal transduction pathway is regulated by multiple genes. Among them, NPR1 (nonexpressor of PR-1) is a key gene located downstream of the SA signal pathway. NPR1 interacts with the transcription factor TGA to activate the expression of a series of SAR genes, and finally confers improved resistance of plants to diseases 71 . Yang et al. 72 reported that enhanced activation of SA biosynthesis in Arabidopsis thaliana hybrids may contribute to their increased resistance to a biotrophic bacterial pathogen. In this study, NPR1 was down-regulated in both R-3dpi vs. R-mock and S-3dpi vs. S-mock comparisons, but PR-1 was up-regulated only in the S-3dpi vs. S-mock ( Supplementary Fig. S3). Moon et al. 73 suggested that overexpression of OsTGA2 can improve the resistance of rice against bacterial leaf blight disease by directly regulating the expression of defense-related genes. Liu et al. 74 found that overexpression of NtPR1a contributed to increasing resistance to R. solanacearum in tobacco Yunyan 87 via activating the defense-related genes. Aux/IAAs (Auxin/indole-3-acetic acid proteins) play important roles in auxin signaling pathways, Fan et al. 75 identified some MeAux/IAAs from Manihot esculenta as novel members in plant disease resistance through gene profiling and functional analysis, and these observations provide important information for further utilization of MeAux/IAAs. Similarly, AUX1 was found down-regulated after inoculation with pathogen ( Supplementary Fig. S3). The plant hormone abscisic acid (ABA) was found in various parts of the plant. ABA can convert the adversity stimulus in the environment into internal signals of plant cells, and further, make plants to produce physiological reactions and resistance towards the invasion of pathogens 76,77 . Earlier studies have shown that the protein phospholipase 2 (PP2C) and sucrose non-glycolytic protein kinase 2 (SnRK) genes in plants can eliminate the inhibitory function of phosphate groups when the content of ABA is low, and then the production and release of ABA will be increased. It further promotes the formation of the ABA-PYR/PYL-PP2C complex, which in turn activates SnRK2. SnRK2 further activates the expression of downstream components, stress-related transcription factors, secondary messengers, and related genes, and regulates stress-responsive gene expression, ultimately protecting plants from adverse stresses [78][79][80][81] . In our study, PP2C and SnRK2 were up-regulated in both R-3dpi vs. R-mock and S-3dpi vs. S-mock comparisons, but PYR/PYL was down-and up-regulated in R-3dpi vs. R-mock and S-3dpi vs. S-mock, respectively ( Supplementary Fig. S3). In summary, several genes in the SA, auxin, ABA, and ET signaling pathways are involved in the defense response against R. solanacearum infection in tomato plants.
In the process of long-term interaction with pathogenic bacteria, plants have formed a set of the natural immune system, which includes two levels, namely PAMP-triggered immunity (PTI) and effector-triggered immunity (ETI). PTI promotes the recognition of PAMPs through cell surface receptors, such as bacterial flagellin and fungal chitin, thereby inducing host plants to produce a series of defense responses, including the formation of phytoalexins, the expression of disease-related proteins, and the deposition of the corpus callosum. ETI means that plants recognize the avirulence gene (Avr) of pathogens through the resistance gene (R), which triggers a series of specific defense responses [82][83][84] . In this study, Rd19 was up-regulated and CDPK, FLS2, and CERK were down-regulated in both R-3dpi vs. R-mock and S-3dpi vs. S-mock ( Supplementary Fig. S2). Shi et al. 85 reported that silencing of TaCP1 (an RD19-like cysteine protease) reduced wheat resistance to Puccinia striiformis f. sp. tritici. Fu et al. 86 showed that rice OsCPK10 is an important regulatory factor in plant immune response, which may regulate disease resistance by activating SA-and JA-dependent defense responses.
In addition, the study also found that the MAPK signal transduction pathway-related genes that mediate plants in response to pathogen infection are up-/down-regulated, respectively, including genes of MAPK and PR-1 (Supplementary Fig. S4). MAPK is a type of serine/threonine-protein kinase, which is ubiquitous in plants. It can receive extracellular or intracellular signals, and further integrate and amplify the signals, ultimately causing cell physiological responses. PR-1 is a pathogenesis-related protein (PRP) that is closely related to disease resistance in plants. This study has also been shown that various stimuli such as pathogen infection, mechanical damage, or pathogen elicitor can activate MAPK, and the MAPK system can further activate related transcription factors. The domains of transcription factors can combine with the W box (T) TGACC (A/T) of PR-1 and can quickly activate the early defense response signal of plants 87 . In this study, genes that displayed up-/downregulated expression indicate that the MAPK signal transduction pathway is involved in the response of tomato to R. solanacearum.

Conclusion
In this study, we analyzed the responses of BW-resistant and BW-susceptible tomato lines to R. solanacearum using RNA-seq. Gene expression changes following R. solanacearum-inoculation were identified, the expression levels and types of DEGs were investigated, and GO and KEGG enrichment analyses were performed to annotate Scientific Reports | (2022) 12:22137 | https://doi.org/10.1038/s41598-022-26693-y www.nature.com/scientificreports/ the function of DEGs. Our results suggest complex and different responses of the two tomato inbred lines to R. solanacearum infection. This study provides an overall representation of the regulatory gene network for both resistant and susceptible tomato inbred lines responding to R. solanacearum infection. These R. solanacearumresponsive genes could serve as important candidates for future functional research and may be helpful in developing an effective method to resist this significant disease.

Materials and methods
Plant materials and growth conditions. Two S. lycopersicum inbred lines (selected and bred in our laboratory), Hm 2-2 (R, resistant to BW) and BY 1-2 (S, susceptible to BW), were used in this study. No approvals were required for this study, and we complied with all relevant regulations during the experiments. The collection, preservation, and use of plant materials involved in this study complied with relevant institutional, national, and international guidelines and legislation. S. lycopersicum seedlings (preserved in our laboratory) were planted in plastic trays (54 × 28 × 5 cm) filled with substrate (peat: perlite = 3:1, v/v).They were placed in an artificial climate chamber with a temperature of 28 ~ 30/15 ~ 17 °C (day/night), a photoperiod of 14/10 h (day/ night), and relative humidity of 50% at Yichun University, Yichun, China. Three weeks later, 100 seedlings at the three-leaf stage were transplanted into black plastic pots (130 mm height × 150 mm diameter) filled with substrate (peat: perlite = 3:1, v/v), and plants were set 10-13 cm from each other. These plants were maintained in an artificial climate chamber with growth conditions as described above. RNA-seq data analysis. Raw reads obtained from the sequencer contain adapters or low-quality bases, which would affect the following assembly and analysis. Thus, to get high-quality clean reads, raw reads were filtered by fastp v0.18.0 89 . Short reads alignment tool Bowtie2 v2.2.8 was used for mapping reads to the ribosome RNA (rRNA) database 90 . The rRNA-mapped reads were removed. The remaining clean reads were further used in assembly and gene abundance calculation. An index of the reference genome was built, and paired-end clean reads were mapped to the S. lycopersicum reference genome (ITAG3.2) using HISAT v2.2.4 with "-rnastrandness RF", while other parameters were set as default 91 . The mapped reads of each sample were assembled by using a reference-based approach by StringTie v1.3.1 92,93 . For each transcription region, an FPKM (fragment per kilobase of transcript per million mapped reads) value was calculated to quantify its expression abundance and variations, using StringTie v1.3.1 software.

Analysis of DEGs.
DEGs were detected using the DESeq2 package 94 . The analysis mainly consisted of three steps: (1) normalization of the read count; (2) calculation of the probability of hypothesis testing (p-value) according to the model; (3) conduction of multiple hypothesis testing to obtain the false discovery rate (FDR) value. The genes with a p-value ≤ 0.05 and |log2FC|≥ 1 were considered DEGs.
GO enrichment analysis. Gene ontology (GO) is an international standardized gene functional classification system that offers a dynamic-updated controlled vocabulary and a strictly defined concept to a comprehensive description of gene properties and their products in any organism 95  Pathway enrichment analysis. Genes usually interact with each other to play roles in certain biological processes. Pathway-based analysis helps to further understand gene biological functions. KEGG (Kyoto Encyclopedia of Genes and Genomes) is the major public pathway-related database 96,97 . The KEGG pathway was taken as the unit and a hypergeometric test was used to identify pathways that were significantly enriched among the DEGs compared to the entire genome background. Through significant enrichment of pathways, the most important biochemical metabolic pathways and signal transduction pathways involving DEGs can be determined.

Single-nucleotide polymorphism (SNP) analysis. Variation detection based on transcriptome
sequencing mainly includes SNPs and insertions/deletions (InDels). SNP variation refers to the DNA sequence polymorphism caused by the change of a single nucleotide at the genome level, while InDel is the variation caused by the insertion or deletion of a nucleotide. The software GATK v3.4-6 was used to detect the variation of SNPs and InDels and finally conduct data statistics 98 .
Validation of RNA-seq data by qRT-PCR. To validate the RNA-seq results, quantitative real-time PCR (qRT-PCR) was performed in this study. RNA was extracted using the HiPure Total RNA Mini Kit (Magen, Guangzhou, China). Later, RNA was reverse transcribed to cDNA using the HiScript II 1st Stand cDNA Synthesis Kit (+ gDNA wiper; Vazyme, Nanjing, China). The qRT-PCR was performed in triplicate for each sample using the 2 × RealStar Green Fast Mixture (with ROX) as per the manufacturer's instruction manual. The qRT-PCR amplification was performed using the StepOne Real-Time PCR Instrument (Applied Biosystems, Thermo Fisher, USA) and corresponding software (Applied Biosystems). The tomato Actin (Solyc03g078400) gene was used as the internal control 99 . The relative expression levels of the genes from three biological replicates were calculated using the 2 −△△Ct method 100 . All primers for qRT-PCR are listed in Supplementary Table S1.

Vector construction and virus-induced gene silencing (VIGS) transformation.
The VIGS vectors pTRV1 and pTRV2 101 were stored in our laboratory (Fig. 9). A 199-bp WRKY75 (initially characterized by López-Galiano et al. 102 and Chen et al. 103 ) fragment was amplified from the stem tissue of resistant tomato plants (Hm 2-2) with the specific primers V-F and V-R (listed in Supplementary Table S2) by PCR. The WRKY75 fragment and the vector pTRV2 were then digested with EcoR I and BamH I (Takara, China), and the WRKY75 fragment was ligated into the vector using T4 ligase (Takara, China) to create the TPV::WRKY75 construct; pTRV2 with empty was created to TRV::empty and was used as control. The resulting vector was introduced into Agrobacterium tumefaciens GV3101 (WEIDI, China). The VIGS experiments were carried out as described previously 104 . Three-week-old newly emerged leaves of Hm 2-2 tomato plants were transformed with Agrobacterium containing TPV::WRKY75 (V(pTRV1):V(pTRV2::WRKY75) = 1:1) and TRV::empty (V(pTRV1):V(pTRV 2::empty) = 1:1) vectors, respectively, using a 1-mL syringe. Each experiment included three biological replicates. One week later, the plants were inoculated with R. solanacearum by using the root-cutting and root-grafting method. Disease scoring was performed according to that previously described by Lacombe et al 105 .
Statistical analysis. Data were expressed as the mean ± standard deviation (SD). Statistical analyses were performed by SPSS (Statistical Product and Service Solutions) and Excel 2007. Differences were considered statistically significant at p < 0.05 and p < 0.01.